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Abstract 

The analysis of scattering from complex objects using surface integral 
equations is a challenging problem. Its resolution has wide ranging 
applications- from crack propagation to diagnostic medicine. The two 
ingredients of any integral equation methodology is the representation 
of the domain and the design of approximation spaces to represent 
physical quantities on the domain. The order of convergence depends 
on both the surface and geometry representation. For instance, most 
surface models are restricted to piecewise fiat or second order tessella- 
tions. Similarly, the most commonly known basis spaces for acoustics 
are piecewise constant functions. What is desirable is a framework that 
permits adaptivity (of size and order) in both geometry and function 
representations. Unlike volumetric, differential equation solvers, such 
as the finite element method, developing an hpadaptive framework for 
surface integral equations is very difficult. This papers proposes a res- 
olution to this problem by developing a novel framework that relies on 
reconstruction of the surface using locally smooth parameterizations, 
and defining partition of unity functions and higher order basis spaces 
on overlapping domains. This permits easy refinement of both the ge- 
ometry and function representation. This capabilities of the proposed 
framework are shown via a number of numerical examples. 

PACS numbers: 43.20.Fn, 43.58.Ta, 43.28. Js 



I. INTRODUCTION 

The analysis of scattering from complex objects has wide spread applicability, from 
crack propagation^l^l to non-destructive evaluatiorpHl to imaging and diagnostic medicine^ 
to holography^ to scattering from rough surfaces^, etc. Boundary integral formulations offer 
an efficient modality for the analysis of fields scattered by homogeneous objects as (i) they 
can be formulated only in terms of surface integral equations and (ii) radiation boundary 
conditions are explicitly included in the Green's function. Despite their advantages, their 
formulation is more difficult than that of their differential equation counterparts, and as a 
result this method has seen sporadic development in the paslP^'^, and a more concerted 
effort recentl} ^^ * ^^ * ^^ * ^^ * ^^ * ^^ * ^^ * ^^ !. The recent development of fast solvers that ameliorate the 
CPU and memory complexity of surface integral equation based solvers, i.e., reduce the 
scaling from 0{N^) to 0{Ns log2 Ng) where Ns is the number of spatial degrees of freedom, 
has made these techniques more appealing^^Ei] However, when compared to their differential 
equation based counterparts, the analysis here has been more or less restricted to simple 
basis functions (piecewise constant) and linear tessellations of the geometry. This is not to 
say that there is not a need for higher order function and geometric representatiorP^E^E^EH. 
Indeed,'^ makes an eloquent case for the development and use of such methods for scalar 
finite element problems. In this paper, our objective is to develop a fiexible framework such 
that both the surface and function representations lend themselves to adaptivity in terms of 
patch size (h—), surface order (g—) and polynomial basis order (p— ). In what follows, we 
shall review extant literature and motivate the need for such a method. 

Constructing an underlying mesh/tessellation is, perhaps, the most understated task in 
any computational analysis. Commercial meshing software exists that provides higher order 
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tessellations, upto second order representation of surfaces. These meshes, while adequate for 
many applications, present challenges when one requires /i— , p— or (?— adaptivity to ensure 
convergence of the solution or flexibility (mesh adaptation) to model crack propagation or 
deformation or the ability to handle non-conformal meshes (meshes that contain hanging 
nodes). Note, that the construction of a framework for h—, p— and (7— adaptivity for volu- 
metric tessellations, such as those used in finite element methods is comparatively easy. It 
is apparent that a possible way to accomplish these goals is to define the topology of the 
scatterer using point clouds, and then use this to create local surface descriptions. Using 
point clouds to define a mesh is not new; it has been extensively studied with regard to effi- 
ciency in two dimensions and 0{N log N) algorithms exisl)^^* ^*^^ * ^^ * ^° l In three dimensions, 
creating a surface description using a set of point clouds is a highly overdetermined problem 
and therefore, considerably harder. A very successful algorithm to create tessellation from 
a set of points is the ball pivoting algorithm (BPA)0^. However, while the BPA has some 
deficiencies, the fundamental problem with this approach is that it creates a mesh. Once 
this mesh is created, it suffers from the same lack of adaptivity and flexibility alluded to 
earlier. Our approach to solving this problem is to develop a methodology that relies on 
defining piecewise overlapping domains. A local polynomial representation of these domains 
is obtained by ensuring that this representation matches the surface at a dense set of points 
specified by the point cloud. Given a local representation in terms of analytic functions, 
mechanisms are then developed to divide or merge patches (/;,— adaptivity), or change the lo- 
cal polynomial order of the surface ((7— adaptivity), or both. Once a surface parametrization 
is obtained, the next step is the construction of function approximations on these surfaces. 
To this end, the two most desirable attributes of the surface parametrization are that: 

1. It be local in nature. In particular, the parametrization should lend itself to the 
definition of local approximation spaces 

2. It permits easy definition of surface derivatives and differential forms. All integral 
equation formulations require the construction of integrals on the surface and some 



formulations require multiple surface derivatives of functions. Thus, it is desirable 
that the surface parametrization allow for the construction of surface derivatives and 
Jacobians in closed-form. 

Our surface parametrization scheme draws inspiration from algorithms that exist in com- 
puter graphics^^'^^'^ and computational physics^ literature. However, to our knowledge, 
the algorithm presented here is the first that satisfies both the properties described above. 

The other critical component of scattering analysis is the approximation space. As 
alluded to in^^, it is eminently desirable to use higher order functions for approximation, 
or better still, use functions that are based on the known local physics/heuristics. As is 
well known, these methods produce higher rates of convergence^^'^. Further, these features 
would reduce computational cost without the detrimentally affecting accuracy. However, 
it is challenging to mix different orders of basis functions, for instance, when one desires 
p— adaptivity or the use physically relevant basis sets in each region. Incorporating such 
flexibility into classical solvers is very difficult as one has to take steps to ensure continuity of 
the physical quantity being represented. In the finite element community, the need/desire to 
enrich the approximation space, as well as ensure continuity gave rise to the generalized finite 
element methocpSEHSlST]^ ^^^ j^g variation d^^ * ^^ * ^° l The basis functions developed within this 
framework are continuous across domains and, as a result, do not need additional constraints 
to ensure continuity. The authors have recently developed methods that extend this idea 
to surface integral equations as applied to electromagneticg^^ES^ an^i have demonstrated 
convergence, well conditioned properties as well as application analysis of scattering from 
the range of targets. Unfortunately, this method still relies on an underlying tessellation. 
While Nystrom based schemes that use a point cloud and collocation, instead of a mesh 
and basis functions, are available^^'^, these methods have been shown to be equivalent to 
higher order basis function schemes on a standard tessellatiorpSESl^ ^g g^ result, both of these 
methods suffer from the drawbacks alluded to in the earlier paragraph. 

The goal of this paper is to introduce a very general and flexible framework for surface 
integral equation based analysis of surface scattering. We will consider acoustic scattering 



from sound-hard objects as the target problem in this work. The method, called the General- 
ized Method of Moments (GMM), aims to introduce this range of functionality. Specifically, 
in this paper we will 

1. present the GMM computational framework that will 

(a) introduce a framework to develop local analytical surface representations on over- 
lapping domains starting from either a tesselated object or a point cloud 

(b) develop the mechanisms necessary for either merging or partitioning subdomains 

(c) develop the mechanisms for locally increasing/decreasing the order of represen- 
tation of each subdomain 

(d) define basis functions with a partition of unity framework that are defined on 
these overlapping domains 

(e) describe accurate evaluation of integrals and the Galerkin solution process 

2. and present results that demonstrate 

(a) /i— , p—, and (7— convergence of surfaces using overlapping subdomians 

(b) convergence of function representation 

(c) convergence of scattering cross-section from canonical geometries and scattering 
cross-sections from topologically different objects 

(d) the ability of the method to mix different basis functions (polynomial or non- 
polynomial) in different regions or basis function adaptivity 

(e) h,p—, hp—, and gf— adaptivity (in both surface and basis functions) 

Note, while it is not a direct focus of this work, the GMM framework introduced here 
can be easily accelerated using the fast multipole methocP^'^ to permit the analysis of very 
large objects. The framework constructed here also permits easy integration with the fast 
multipole method The rest of this paper proceeds as follows: In the next section, we will 



formally state the scattering problem. In Section III , we will describe the construction of the 



locally smooth surface parametrizations using either an underlying mesh or a point cloud. 



Section IV details the construction of the basis functions on these surface parametrizations. 



The specifics of construction of the matrix elements will be elucidated in Section IV] and Sec- 



tion VI will present several results that demonstrate the surface reconstruction, validate the 



basis function framework and showcase the advantages of the proposed technique. Finally, 



Section VII will provide some concluding remarks. 



II. PROBLEM STATEMENT 

Let D~ denote a rigid scatterer in a homogeneous medium bounded by Vt with a unique, 
outward pointing normal n(r)Vr G i7. Consider a velocity field incident on this scatterer 
denoted by v*(r). This generates a scattered velocity field given by v*(r) and we define the 
total velocity as v*(r) = v'(r) + v''(r). These fields can be represented by an equivalent 
potentials 0^(r), C, € {z,s,t}, where v''r) = V0^(r). Further, the corresponding pressure 
fields are given by p'^ij) = ~jujpQ(f)^{v) where po is the density of the ambient medium. The 
total potential 0*(r) = 0*(r) + 0*(r) satisfies the Helmholtz equation and boundary condition 
given by 

VV*(r) + A;2</)*(r) = V re TZ^ / D' 

(1) 
n(r) ■ V(/)*(r) = V r eVt. 

The Kirchoff- Helmholtz integral theorem relates the scattered potential 0^(r) to the total 

potential as 

<P'{r) = I rfr0*(r')n'(r') ■ V'^(r, r'), (2) 

Jvt 

where g{v^ r') = exp(— jA;|r — r'|)/47r|r — r'| and k is the wave number of the incident field. 
Imposing the condition that the total pressure p*(r) = p*(r) +p*(r) = on the surface Vt 
provides an integral equation for the total potential, 0*(r), given by 

<t>\v) = \<p\v) - I rfr>*(r')n'(r') ■ V'^(r, r'). (3) 



Further, by imposing that the normal component of the velocity goes to zero on the 
surface of the scatterer, i.e. n(r) ■ v* = 0, V r G ^2, we obtain the normal derivative of the 
above integral equation. 

n(r) • V0^(r) = / rfr>*(r')n(r) • Vn (r') ■ V'^(r, r'). (4) 

Jn 

We define two integral operators /C and T as 

/C o [0(r)] = ^0(r) - f rfr>(r')n (r') ■ V'^(r, r') (5a) 

2 Jn 

and 

r o [(/)(r)] = f (irV(r')n(r) ■ Vn'(r') • V'^(r, r') (5b) 

Jn 

The two integral equations in (Is]) and ^ can be combined using a parameter a as follows, 
in a formulation that guarantees uniqueness in the solution (/)*(r)'ni 

a(j)'{r) + (1 - «)n(r) ■ V</>^(r) = aJC o [0*(r)] + (1 - a)T o [0*(r)], (6) 

where a G (0, 1). Solution of equation ^ by the method of moments proceeds by represent- 
ing the unknown potential 0*(r) in a set of spatial basis functions, i.e. 0*(r) = X]n'^"*^"('^)' 
where a„ are unknown coefficients. Substituting this representation into (|6| and using 
Galerkin testing results in a matrix system of the form 

Za = f, (7) 

where 

Z = [^d= / t/r0,(r)A'o[0^.(r)], (8) 

and X = a/C + (1 — a)T, a = [oj] and 

f =[/,]= f dr<jy,{v)<jy'{v). (9) 

Typical method of moments solutions employ polynomial basis functions defined on a simpli- 
cial tessellation of the geometry ^2. These basis spaces rely on mapped polynomial functions 



defined on each simplex. In keeping with the goals stated in Section [T| our approach to 
solving this problem is to (i) develop an overlapping local surface parameterization each of 
which will be the domain of the support of the basis function and (ii) define basis functions 
on these domains. Insofar as the latter is concerned, it builds upon the framework developed 
jjj 47 | 48 | 4i |42] ^^^ piecewise fiat domains. 

III. CONSTRUCTION OF LOCALLY SMOOTH SURFACE 
PARAMETRIZATIONS 

Next, we prescribe the construction of domains of support for basis functions; these 
domains overlap and form an open cover of the scatterer Q. For the rest of this section, we 
will assume that we have a point cloud, a set of normals to the surface at these points and a 
connectivity map that identifies nearest neighbors for each point. Note, that algorithms such 
as ball pivoting^ may be used to obtain such a connectivity map in linear time. Algorithm 
[l] presents a sequence of tasks in order to construct these patches. The rest of this section 
will elucidate each of the steps presented therein. 

Algorithm 1 (Color online) Outline of patch construction 
1: Subdivide initial primitives into overlapping neighborhoods {fi,} 

2: for each neighborhood Qi do 

3: Project each neighborhood onto a plane Fj 

4: Construct a local coordinate system {u, v, w) using the projection plane and its normal 

5: Construct GMM patches Aj as a least squares, polynomial approximation to Qi 

6: Merge or split these patches if necessary 

7: end for 

A. Construction of GMM domains 

We begin by partitioning the domain Q into neighborhoods Qi that overlap and 
completely cover the domain. To this end, assume that the domain Q is described by a set 



of nodes Nl = Ufl^jA/i}, a connectivity map consisting of primitives A^v = U^^]^{A„}, and 
finally a unique set of normals fij at these points. Each primitive is defined by a collection 
of nodes A„ = {A/'„j}"'~™" C Ml- In the case of a standard, flat, triangulation, this will 
reduce to m„ = 3 Vn, i.e., all the primitives are triangles. To define locally smooth GMM 
patches, we first start from a collection of primitives that share a node Mi. This collection 
will be denoted by fij, called the GMM neighborhood. A set of neighborhoods constructed 
from a point cloud is described in Figure [l} To construct a locally smooth approximation 
to VL starting from these neighborhoods we first define a neighborhood normal, a projection 
plane and a notion of permissibility for each neighborhood Vti as follows: 
Definition 1: Permissible neighborhoods 

Given a neighborhood fij, centered around a point Mi = r^ and a parameter e, the 
average normal for the neighborhood Qi is defined as 

^. = ]^E^EftM. (10) 

n=l fc=l 

where Ni is the number of primitives A„ connected to A/i, r^ are points chosen on each of 
the member primitives such that r^ C fij fl A^, and \\Yk — i"«||2 ^ ^• 

Further, a projection plane for neighborhood Vti is defined as the plane passing through 
Yi and normal to fij^e. Let Fj be the projection of Vti on this plane and denote the projection 
of a point r G f2j to the plane Fj by r'. The neighborhood Vti is permissible if we can find 
some e such that, Vr G fij and for r 7^ r^, 

> 0. (11) 



(r - Yi) X ni,e 



In other words, a permissible neighborhood is one for which we can find a projection 
plane such that the entire neighborhood lies on one side of the plane and has a unique 
projection on the plane. Figure |2] shows the construction of the neighborhood normal and 
projection plane for a permissible neighborhood. For each permissible neighborhood, we 
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FIG. 1: (Color online) GMM neighborhoods ({^j}) constructed by partitioning a point 
cloud shown as shaded region. The neighborhoods are constructed as a set of nodes 

connected to a node. 





FIG. 2: (Color online) Construction of a projection plane (Fj) from a GMM neighborhood. 



can define a local coordinate system containing the projection plane and the neighborhood 
normal, as follows: 



Definition 2: Local coordinate system 

For a permissible neighborhood Qi, choose a point r'^ such that |r'^ — r^| > 0, on 
Fj and define the following local co-ordinate system {u, v, w}j and corresponding pro- 
jections u(r), v(r),w(r) for any point r E Qi as w = hi, u = (r'^ — r'-)/|r^ — r^l, 
V = (ii X u)/|n X "ii|, u(r) = r' ■ u, v(r) = r' ■ v and w(r) = r' ■ w. 

Finally, the above definitions can be used to generate a polynomial map whose domain 
is the projection Fj, described by the local coordinates (u, v) and whose range is a smooth 
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FIG. 3: (Color online) Construction of a local coordinate system on the projection plane. 



surface. This mapping will become the "generator" for the locally smooth surface and is 
called the GMM surface map. 



Definition 3: GMM surface map 



Given a permissible neighborhood Qi and a corresponding coordinate system {u, v, w}, 
we can define a polynomial "Pf (u, v) in two variables (u,v) complete to order g by its coeffi- 
cient vector Cf = [cq, . . . C(^g+i){g+2)/2\ ■ The polynomial Vi{u, v) (and corresponding C^), that 
minimizes the norm min ||Pf (u(r), v(r)) — w(r)|L can be used to define a transformation 
£f from Qi to Aj, given by 

£f (r) : a ^ Ai = uu + vv + Pf (u, v)w. (12) 

Aj forms an order-p smooth, least-squares approximation to fli. This transformation will 
be called the GMM surface map. The patch Aj will be called a GMM patch of order g. 
In typical implementations, the "user" either chooses a desired order of the patch gr or an 
error criteria Er- Correspondingly the error Eg or order g^ is determined by the minimization 
procedure. The error is computed at a random selection of points on the patch. 
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FIG. 4: (Color online) Construction of a locally smooth parameterization starting from a 

point cloud. 

B. Merging and splitting of patches 

Next, to achieve complete h— and gf— adaptivity of the surface parametrization, it is 
necessary to construct a scheme to merge or subdivide patches as necessary. In this section we 
will detail "adaptive" algorithms to merge or split patches based on a "sharpness" criterion. 

1. Merging patches 

Let Vti and VLj be two GMM neighborhoods with average normals fij ^ and nj_£ respec- 
tively. The neighborhoods are said to fail the smoothness criterion if fij ^ ■ fij ^ > Smi where 
Em is a user-determined smoothness threshold. The two smooth neighborhoods are merged 
such that fifc = Vti IJfij, and a new smooth parametrization A^ is constructed from Vt^. 

2. Splitting patches 

Let Vti be a GMM neighborhood with average normal hj^£. For each point A4 = i"fc in the 
neighborhood, the neighborhood is said to fail the sharpness criterion at Mk-, if n(rfc) ■ fij ^ < 
Esi where Eg is another user-determined sharpness criterion. If a neighborhood fails the 
sharpness criterion at A/fc, the point Mk is excluded from the neighborhood Vti and a new 
neighborhood fi^, is constructed using primitives that share A4- Further, two new GMM 
patches Aj and A^ are constructed using (the new ) neighborhoods Vti and Vtk- Note that, 

1. Since the merging and splitting are operations constructed at the neighborhood level, 
these operations can be performed either before or after patch construction. This is 
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important if the aim is to obtain /;,— adapt ivity. 

2. The merging and sphtting processes can be recursively repeated until all neighborhoods 
pass both the smoothness and sharpness criteria. 

3. The smoothness and sharpness criteria need not be global to the problem; in fact, 
depending on the level of control desired, these criteria can even be on a "patch-by- 
patch" basis. 

4. If the aim is to obtain a true /ip— refinement scheme, solutions can be constructed and 
the patches can be merged/split as necessary to refine the solution. 

C. Local derivatives, normals and continuity of functions 

In order to construct functions and surface derivatives on the locally smooth patches, we 
need to construct surface gradient tensors. Given the GMM surface map Cf{r), we denote 
its first metric tensor by 



G. = 



9u 912 

921 922 



duV ■ duT duT ■ dyV 
(9„r ■ a,r a,r ■ d„r 



(13) 



The corresponding surface differential element is denoted by 



dS = ^/(fidudv, 



(14) 



where gi = det{Gj), the determinant of the metric tensor. Each term in the tensor can be 
defined in terms of the polynomial Vf{u,\) as 



(15) 



duV = u + duVfiu, v)w, 
(9„r = V + d^V-{u, v)w. 

Given a scalar function 0(u, v) defined on the projection plane Fj, the surface gradient 
of the function on Aj is given by 



V^0 = 9ndu(l) duY + 5'i2<9«0 d^r + 5'2i<9„0 <9„r + 5'22<9t,0 d^r. 



(16) 



14 



Higher order derivative tensors on the surface can be described in a similar manner. Finally, 
the surface normal at any point on the GMM patch can be defined as 

duvxd^r 

Note, that this normal is continuous across the entire patch up to one less than the order of 
the patch g. 

From the definition of the surface gradient above, it is clear that any function 0(r) 
defined on a GMM patch Aj of order g that supports p derivatives on (u, v) with p < g 
will support at least p surface derivatives on the smooth patch Aj. This result implies that 
defining a function of order p on the smooth GMM patch Aj is equivalent to defining a 
corresponding function on the projection plane Fj. This provides an important tool for 
defining GMM basis functions as described below. 



IV. DEFINITION OF GMM BASIS FUNCTIONS 

The next step is the development of basis functions in each of the above patches. Con- 
sistent with the central theme of the GMM framework, we develop a scheme that permits 
different orders of polynomials or different functions to be defined on adjacent patches. It 
has been shown, for integral equation based solvera^^ESl^ ^j^g^^^ ^j-^jg ^g^y^ i^g achieved using a 
product of two functions; (i) a partition of unity (PU) function that provides continuity of 
the order of this function across overlapping patches and (ii) a higher order function that 
determines the quality of approximation within a patch. In what follows, we shall briefly 
discuss each in turn for completeness. Details of development of basis functions can be ob- 
tained fromP^ESl^ with sufficient modification so as to include the general nature of the local 
surface description. 
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FIG. 5: (Color online) Definition of a GMM patch and partition of unity 

A. Definition of partition of unity functions 

Consider a GMM patch Aq. Let Aj overlap with Ni other patches {Ajjj^^. Then a 
partition of unity function is defined on Aj as a function ^j(r) that satisfies the following 
properties 



1. ^i(r) = V r ^ A, 

2. ^P,{r) + E7=o^.(r) = 1 V r G J (a. {A,}f=i 



In practice, to define a partition of unity function on Aj, we construct a function Aj(u, v) 
which is 1 at the patch center and at the edge of Fj, the projection of Aj. The partition 
of unity is then defined as 

where the index k runs through all the patches Qk that overlap with fij. It can be verified 
that this definition ensures that the partition of unity goes to at the ends of the patches and 
adds up to 1 everywhere on Fj. Correspondingly it satisfies these properties on Aj. Higher 
order PU functions can be defined in a similar manner if necessary. For illustration. Figure [5] 
shows two one-dimensional patches and a partition of unity defined on these patches. Figure 
[6] shows the construction of Aj j for a fiat patch. 
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FIG. 6: (Color online) Definition of a pyramid function for partition of unity - A 
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B. Definition of continuous approximation functions 

The next step is to define functions that provide higher order approximation of the 
unknown field in the patch. As before, we start by defining the function on Fj. Any function 
/(u, v) can be now mapped directly to /(r) on Aj. Note, the domain of the approximation 
function does not need to be identical to the projection of the patch, Fj. This is possible as 
functions defined on these patches are eventually multiplied by a PU function that goes to 
zero at patch boundaries. 

One possible choice of approximation functions can be described using Legendre polyno- 
mials of the form /^"(r) G {Pp^{'u)Pp^{v)} where Pg denotes a Legendre polynomial of order 
q and p^ + p^ < m and 

~, ^ . uir) 

ulr) 



maXreA^ u{r) 



v{r) = ^i^. (19) 

maXreA, ^^(r) 

Once approximation functions are thus defined, the GMM basis functions are simply prod- 
ucts of the approximation function with the partition of unity. That is, 

<Pi{r)espan„,{^i{r)ur{r)} (20) 

Once the basis functions are defined, the next step is the evaluation of the integrals to 
construct the matrix elements in [Zij]. This will be detailed in the next section. 
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V. EVALUATION OF MATRIX ELEMENTS 

The evaluation of the matrix elements in [Zij] involves integrals of the following two 
forms. 

/ rfr0,(r) / dr'<P,{r')h'{r')-Vg{r,r') (21) 

JAi JAj 

[ dv(t>i{v) [ rfrV,(r')n(r) ■ Vn'(r') ■ V'^(r,r'), . (22) 

JAi JAj 

The integrals need to be evaluated on patches Aj and A,-. Using the surface differential 



element defined in (14), we can map the integral of a function B(r, r') on a patch Aj to an 



integral of the function B(u,v,u', v') on the projections Tj and Fj as 

dv dv'(p{r) = / ^/oidudv I y/gjdu'dv' Q{u,\^,u',V). (23) 

JAi Ja, JVi JVi 



The evaluation of the integrals in (21) and (22) are performed using the transformation 
in (23) and Gaussian quadrature when the patches are well separated from each other. It 
is observed the Gaussian quadrature rules converge to sufficient accuracy when the centers 
of the patches are separated hj d > 0.15A, where A is the wavelength of the incident field. 
When the patches are closer to each other, the integrals need to be handled more carefully. 
We separate the "near" evaluations into two cases. 

1. Aj and Aj are closer than 0.15 A but do not overlap: In this case, the integrals are near 
singular, but can be evaluated using the techniques described iiJUEMSl. 

2. Aj and A^ overlap : In this case, we split the projections Fj and Tj into an overlapping 
section F° and two non overlapping sections Fj/F° and Tj/T". Any integral of the 



form (23) above can be then re- written as follows 



dudv / du'dv' = / dudv / du'dv' + / dudv I du'dv' 
Ti Jr, ir,;/r° ir,/r° ir,/r° ir° , , 

+ / dudv / du'dv' + / dudv / du'dv'. 
JTj/T° jv° jt° jv° 

The preceding equation contains three double integrals that are near singular and one 
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over A" that is either singular (for (21)) or hypersingular (for (22)). The near singular 
integrals are handled as in case 1 above. To evaluate the singular integrals, we make 
the assumption that the overlapping portion is locally flat. This implies that ^ = 1. 



In this case, it can be shown that the integral in (21) reduces to 0. The integral of 



equation (22) on flat patches can be performed by transforming the surface integral 



into a line integral as described irt^^. 

VI. RESULTS 

In this section, we present a series of results that demonstrate the features of the GMM 
scheme presented here. The results can be broadly categorized into two: (i) geometry and 
function representation, and (ii) using this representation in Galerkin framework to solve 
integral equations. To start out, we will define error metrics that one can use to define 
accuracy of reconstruction of a surface. This is then followed by results that demonstrate 
the following: (i) convergence of surface reconstruction of analytically describable surfaces 
represented using point clouds; (ii) adaptivity in space and order of surface representation; 
(iii) convergence of function representation; (iv) the ability of the method to use different 
basis functions; (iv) ease of h—, p— , and /ip— adaptivity, and (vi) the capability of analyzing 
geometries described using only point clouds (thus obviating the issues with non-conformal 
tesselations). In all comparisons, we will use scattering cross section (SCS) data that is 
obtained with from analytical results or from a over discretized piecewise constant (function 
and geometry) method of moments solver. 

A. Geometry representation and adaptivity 

In this section, we deal exclusively with various aspects of representing the geometry 
using locally smooth functions, as well as h— (space) and g— (order) convergence of these 
patches. To aid in defining these operations, we define error metrics that will guide adap- 
tivity. 
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1. Error definitions and convergence metrics 

To begin, we present two error metrics that can be used to determine the quahty of 
representation of a surface. These metrics are suitable for defining functions j(r) G H^/'^{VL) 
on the patch A = |J^ Aj and are: 
Definition 4: Surface Approximation Error 

Given a surface approximation A to a true surface Q, the surface approximation error is 
defined as 



^5^||n,(r)nn(r)-nA,(r)||2;£i/2 = ^ fj^ ^ dr n,(r)t(r) - ^ dr t(r) j + £v 



'^ N 

(25) 

where t(r) is any test function and nn(r) and nAi(r) are surface normals to fl and A at 

r G fi and r G Aj respectively; nj(r) is defined by 

1 Vr G QL 
n.(r) = { '' (26) 

else 

Figure [7] demonstrates the convergence of this error on the surface of a sphere of radius Im. 
In order to study convergence, a locally smooth parametrization is constructed starting from 
a two different point clouds. The first, ptri = 1, corresponds to a distribution of points such 
that the average separation distance between nearest neighbors is approximately 0.1m, and 
second, ptri = 2, corresponds to approximately 0.2m. The errors ey and £1/2 are examined 
as a function of the polynomial order of the patch. As is clear from the image, the error 
converges very rapidly with the order of the local parametrization. 



2. Complex geometry representation 

Next, we present results that demonstrate the construction of a surface representation 
directly from a point cloud. It will become apparent that the techniques presented can be 
easily modified to create a smooth surface representation when starting from an underlying 
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FIG. 7: (Color online) Convergence of surface error metrics with surface patch order. 



mesh. To illustrate such a construction, two candidate structures shown are a gyroid and 
an icosahedral geometry enclosing a sphere. 

First, Figure Isla) shows the surface rendering of a gyroid, mathematically described 
by the equation cos(x) sin(|/) + cos{y) sin(x) + cos{z) sin(x) = 1. The surface is complex, 
but as it is analytically known, obtaining a point cloud and corresponding normals at each 
point is relatively simple. Figure [Stb) shows a point cloud constructed from the gyroid 
surface description using — vr < {x, y, z} < n and Npts = 4888 points. Figure Islc) shows 
the results of standard meshing algorithm (ball reconstructiorP^'^ used to create a mesh 
from the underlying point cloud. As is clear from the inset, the resulting triangulated 
mesh has several discrepancies, making it impossible to use in integral equation solvers. 
Furthermore, even if one were to spend sufficient time in cleaning up mesh, it will result 
in systems with high condition numbers as the surface discretization is highly non-uniform. 
Figure IsVd) shows the surface parametrization algorithm that is described in this paper 
applied to the gyroid surface. To construct the parametrization, a primitive is constructed 
at each point in the point cloud and smooth patches are constructed to an relative error 
threshold of Er = 10^'^. As is clear from the figure, it is possible to obtain a locally smooth 
parametrization of the surface starting from a simple point cloud. 

Next, Figure |9](a) shows a point cloud description of an icosahedron enclosing a sphere. 
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(a) Surface rendering of a gyroid 



(b) Point cloud for a gyroid 




(c) BPA algorithm applied to the point cloud 



(d) Construction of locally smooth patches 



FIG. 8: (Color online) Construction of locally smooth surface representation for a gyroid, 

starting from a point cloud. 

Tfie surface of the icosahedron can be constructed in closed form as cos{x + (1 + ^y)) + 



V5, 



V5 



V5. 



V5, 



cos(x-(l + Yl/))+cos(|/ + (l + Y-s)) + cos(x-(y + Y-2)) + cos(z + (l + Ya;))+cos(x-(;z + 
^x)) = 2. Figure 9[b) shows a smooth surface parametrization of the surface constructed 
from the point cloud representation. The point cloud is constructed for an icosahedron of 
radius 3.0 and sphere of radius 0.5 using Npt^ = 4328 points. The final surface is constructed 
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(a) Surface rendering of a icosahedron 
sphere 




(b) Surface reconstruction 



FIG. 9: (Color online) Construction of locally smooth surface representation for an 

icosahedron + sphere geometry. 



to maintain an relative error of e,- = 10 ^. 



Finally, Figure 10 demonstrates the reconstruction of two locally smooth GMM patches 
starting from a piecewise continuous triangulation. Two neighborhoods are defined using 6 
triangles each, shown by Qi and the locally smooth patches Aj are constructed as approx- 
imations to Qi. The error in the patches is computed at 300 arbitrarily chosen points on 
each triangle making up the patch (1800 points overall). In both cases, the relative error is 
maintained to a threshold of Er = 10~^. Note, that the figure is rendered with an artificial 
distance between the flat and smooth parametrization for ease of visualization. The three 
examples provided demonstrate the surface parametrization scheme developed in this paper 
and its ability to construct patches starting from point cloud descriptions of complex objects 
and from a standard piecewise triangulation. 
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FIG. 10: (Color online) Locally smooth overlapping patches constructed from piecewise 

smooth neighborhoods. 

3. Automatic adaptivity of geometric order (g—adaptivity) 

Next, we demonstrate the polynomial adaptivity of the GMM patch. The use of the least 
squares minimization, in the construction of the GMM patch, implies that the polynomial 
order of the map {g) can be automatically chosen depending on an error metric, as opposed 
to being set (by a user) a-priori. To test this property, we consider the error in the surface 
normal to the reconstructed surface, ||ni(r) — n(r)||2. The error is computed with respect 
to the original surface normal at each point in the neighborhood from which the patch is 



constructed. Figure 11 shows the error in the surface normal as a function of the order g, for 
various surfaces, of the form x^" + y^° + z^" = c for a constant c and order parameter go. In 
each case, the surface is first approximated using a point cloud description and then GMM 
neighborhoods Qi are constructed from these triangles. The point clouds are constructed by 
varying x and y in the interval — 1 < {x,y} < 1 and computing z using root finding by an 
exhaustive search algorithm. In each case, the algorithm is run until Npts = 625 valid points 
are found. 

A locally smooth patch Aj is then defined for a given order p and the error in the norm 
is computed. The error convergence in shown for three surfaces, a fiat surface (represented 
as (7o = 0), a piece of a spherical surface {go = 2) and a surface with go = 4. For the latter 
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FIG. 11: (Color online) Convergence of error in normal for surface approximations. 

two cases c = 16.0. In each case, the "true" surface normal can be computed using the 
definition of the surface. As can be seen from the figure, the error for each surface reaches 
machine precision once the mapping order crosses a threshold. This provides a naturally 
adaptive mechanism for the choice of surface order. 



4- Merging of patches (h— adaptivity) 

Next, we present results on the merging of patches based on a smoothness criterion. 
Figures [l2|(a) and fl2|b) show two views of a parabolic surface on which a set of trian- 
gular neighborhoods is constructed. The surface is obtained by constructing the parabola 



y' 



4900x and then rotating it around the x axis. The surface is meshed using Nq = 19930 



triangles. The GMM algorithm is then used to (i) construct neighborhoods using the nodes 
provided by the mesh, (ii) merges these neighborhoods using a constraint on the normals as 



described in Section III.B above, and (iii) construct patches from these neighborhoods. Fig- 
ures 



12 c) - 12 f) show the patches constructed starting from the same initial triangulation, 



using three different thresholds on the angle between the normals (e^ = {5°, 10°, 15°, 22°}). 
In each case, for clarity of representation, only patches that have been constructed from 
more than 500 merged primitives are shown. The patches are constructed to maintain a 
maximum order oi gr = 2 . Table [l] shows the final number of patches {Npat), the maximum 
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order of patches (gpat), the maximum error in the patches Eg and number of primitives in the 
largest patch (Nprim) for each case in the figure. The figures clearly demonstrate the ability 
of the GMM scheme to automatically merge neighborhoods and create smooth patches from 
these neighborhoods. 



B. Function representation 



Next, we will demonstrate the ability of the GMM scheme to represent functions on the 



surface parametrization. Figure [13] shows the convergence of the GMM approximation to a 
function defined on a spherical surface. To test the efficacy of the GMM basis functions, we 
define a function of the form /(r) = f{6,(j)). We then construct local surface parametriza- 
tions {Aj} of varying order f^ = 1, 2, 3 to approximate the surface of the sphere (radius l.Om) 
and construct basis functions of various orders p = 1,2,3 on these surfaces. The functions 
are used to approximate /(r) by setting up the system of equations below. 

/(~r) = ^a,,0,(r) (27a) 

i 

and solving the matrix system resulting from 

/ rfr0, (r) V a,,0i(r) ^ / rfr0,(r)/(r) (27b) 

Jn ^ Jn 

The coefficients Oj are used to approximate /(r) and the norm of the error on the surface is 



used as a parameter to test convergence. Figure 13 shows the error for f(9,(f)) = 6 + (f)as 
a function of p and g, when the patches are constructed from a point cloud on the sphere, 
consisting of Npts = 256 points. A patch is constructed around each point in the point cloud. 
As is clear from the figure, the error converges uniformly and logarithmically with both g 
for each basis function order p. 
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(a) Surface rendering of a parabolic reflector antenna 



(b) Side view 




(c) Patch construction using em — 5" 



(d) Patch construction using em — 10° 




(e) Patch construction using em = 15° 



(f) Patch construction using em = 22" 



FIG. 12: (Color online) Merging patches for automatic h— adaptivity for a parabolic 
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surface. Patches are shown for different merging criteria. 




FIG. 13: (Color online) Error convergence for surface functions defined on a sphere. 

C. Application of GMM construct to solving surface lEs 

Thus far, we have presented results showcasing the surface parametrization and rep- 
resentation of functions on the surface. Next, we use these basis functions and surface 
approximations within the Galerkin framework to solve Eqns. ([6]). To validate the accuracy 
and utility of the GMM technique implemented on the locally smooth surfaces, we preform 
a series of numerical experiments. We begin by presenting results that validate the tech- 
nique on canonical (or near-canonical) geometries. The data obtained using the GMM is 
compared against (i) analytical data and (i) a method of moments integral solver that uses 
fiat triangulation and piecewise constant basis functions. Following this, we will present a 
variety of results that demonstrate (i) h—, p— and /ip— convergence of the GMM scheme, 
(ii) the ability of the GMM to mix different orders and classes of basis functions and (iii) its 
ability to handle complex multiply connected geometries, starting from point clouds. Unless 
specified otherwise, the following is criteria is true for all cases examined: 



1. The surface representation is constructed using a point cloud using and a connectivity 
map such that the distance between each pair of neighboring points is approximately 
O.IA, where A is the wavelength of the incident field. This allows for comparison 
against codes constructed on a standard tessellation with average edge length O.IA. 
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2. The smooth-surface approximations are constructed starting from the point cloud 
provided by this discretization. 

3. In each case, the average radius of the smallest circle containing the projection of the 
GMM patch is used as a measure of the size of the patch, and is maintained to O.IA. 

4. Patches are constructed so as to maintain an error threshold of Er = 10^^. 

In each case, the bistatic scattering cross section (SCS) is used as a metric for comparison 
unless otherwise specified. All the cases demonstrated below assume that the test objects 
are sound-hard and are immersed in a homogeneous medium. The speed of sound in the 
ambient medium is assumed to be 343m/ s. 

1. Validation against analytical results 

First, we consider scattering from acoustically hard spheres of radii O.IA, 0.3A and l.OA. 
The incident velocity field has a frequency of 34.3ifz, and propagates along — z direction. 
The GMM discretizations, result in Ngmm = 300, 450, 500 unknowns for each of the spheres 
when using first order Legendre polynomials p. In each case, the bistatic SCS evaluated at 
= is shown in Figure [l4|^a) and demonstrates excellent agreement with analytical data. 



2. p— and g— convergence 

Next, we consider relative error convergence in backscatter from a sphere of radius O.IA 
between an analytical and the GMM results as a function of (i) the polynomial order of 
the basis functions and (ii) the order of local smoothness of the geometry. The incident 
velocity field is a plane wave of frequency 343Hz propagating along — z. We consider the 
convergence of the relative error in backscatter (0 = 0, ^ = 0). The dashed curve in Figure 
[l4Fb) demonstrates p convergence for fixed h = O.IX and g = 2. The corresponding number 
of unknowns is Ngmm = 320,640,960 for p = 0, 1,2, respectively. As is evident from this 
graph, the error decreases exponentially with increase in p. 
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Next, the solid line in Figure [l4[b) shows the convergence in relative backscatter error 
with the order of the geometry g. The initial patch size is maintained at h = O.IA and the 
polynomial order of the basis functions at p = 1, and as a result the number of unknowns 
is constant at Ngmm = 320. As is clear from the figure, the error decreases exponentially 
with geometry order. 

3. Validation against piecewise constant MoM 

Next, we consider two non-canonical geometries - a NASA almond and a conesphere. 



Figure 15 ^a) shows the bistatic SCS (evaluated at = and = vr/2) due to scattering 
from a NASA almond, that fits in a box of size 3.0A x l.OA x O.IA. A 3A3Hz velocity field 
is incident along — z and the almond is discretized using Ngmm = 1700 unknowns. Figure 
[Isj^b) shows the bistatic SCS (computed at = 0) obtained due to a velocity field incident 
along z on a conesphere with cone-height 2.6A and sphere radius 0.5A. The number of 
unknowns used to discretize Nqmm = 1078. The SCS obtained using the GMM is compared 
against those obtained using classical MoM. Excellent agreement is observed between these 
two data sets. 

D. Mixtures of basis functions 

The results presented thus far validate the GMM scheme and demonstrate convergence 
for various parameters {h, p and g). In what follows, we will present results that demonstrate 
the ease with which different orders and types of basis functions can be mixed together in 
the GMM scheme. This capability is thanks in large part to the fact that the basis functions 
have built in continuity, obviating the need for additional constraints. 

First, consider scattering from an ellipsoid of axes l.OA, 0.5A and 0.25A. The ellipsoid is 
discretized using patches of average radius 0.075A, and the geometry order is maintained at 
g = 2 ioT all the patches. Polynomial basis functions of order p = 1 are used in all patches 
except patches within 0.2A of the two ends of the ellipse. In the patches near the end, 
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(a) SCS due to scattering from spheres of varying radii 
(0.2m, 0.4m, l.Om) 



(b) Convergence of error in backscatter 



FIG. 14: (Color online) Validation results for the GMM -comparison against analytical 
results. SCS from a sound-hard sphere is presented as a metric for comparison 




200 150 -100 



(a) SCS due to scattering from a NASA almond (3.0A x(b) SCS due to scattering from a cone-sphere (2.6A x 
l.OA X 0-1 A) 0-5A) 

FIG. 15: (Color online) Validation results for the GMM - comparison against MoM: 
incident field along z for almond, and x for conesphere. The SCS evaluated at = for 

almond and ^ = for cone-sphere. 



radial basis functions, inspired by^^'^ are functions of the form /( 



u,v 



exp — Cj(m^ + v'^), 



where u, v are the local coordinates on the projection plane, as described in|III[ and Cj is an 
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(a) Basis function partition for ellipsoid 



(b) SCS comparison using multiple basis functions 



FIG. 16: (Color online) SCS computed using mixed order basis functions on surface of an 

ellipsoid. 



order measure, maintained at q = 0.5 for this test. Figure 16 a) shows the partitioning of 



basis functions on the ellipse. Figure 16 b) shows the SCS obtained using this scheme with 



mixed basis functions compared against an SCS obtained using polynomial basis functions 
everywhere {p = 1), and one using radial basis functions (q = 0.5) everywhere. The SCS 
is obtained due to a plane wave incident along x and evaluated at 6' = 7r/2. Excellent 
agreement between all these data sets attest to the flexibility of GMM approach. 

E. hp-adaptivity 

Next, we utilize the flexibility of the GMM scheme to study the /ip— convergence of the 
SCS due to scattering from an ogive of size 10m x 2m x 10m. In each of the cases that 
follows, the SCS is obtained due to a plane wave incident along z, of frequency 3A3Hz. The 
bistatic SCS is evaluated at = 0. To obtain a reference, the ogive is discretized using a 
standard triangulation, at h = 0.05A everywhere and the SCS is computed using a classical 
approach using NmoM = 8406 unknowns. The order of surface parameterization used is 
g = 2 in the smooth areas, (7 = 4 near the ends of the ogive (within 0.25A of the end) and 
g = 7 ioT the two patches near the tips. Using these geometry specifications, the following 
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FIG. 17: (Color online) /ip— adaptivity on an ogival surface: Surface currents and SCS 
comparison for mixtures of various orders of basis functions and patch sizes. 



discretizations are used. First, the ogive is discretized using patches of size 0.25A in the 
smooth areas and 0.1 A near the tips (patches within a sphere of 0.2 A near the tips). Basis 
functions of order p = 1 are used in all patches resulting in Nqmm = 1600 unknowns. This 



case is referred to as hp — in Figure [T7j Next, basis functions of order p = 1 are used in 
the smaller patches and p = 2 in the larger patches. This case is referred to in Figure 17 as 
hp—1, and result sin Ngmm = 2156 unknowns. Finally, the tip of ogive is discretized at O.IA, 
the region near the smooth end of the almond (patches within 0.2A of the smooth end) is 
discretized at 0.15A and the central, smooth portion is discretized at 0.25A. Basis functions 
of polynomial order p = 1, p = 2 and p = 3 are used in each of the areas, respectively. This 
case is referred to as hp—2, and results in Ngmm = 2876 unknowns. The agreement of the 



three different sets of SCS data is shown in Figure 17 and demonstrates the ease with which 
hp convergence can be obtained using GMM. 



F. Application to objects described using point clouds 

Finally, GMM is used to compute scattering from a complex structure that is difficult 
to mesh using any standard meshing technique. Figure [TS] demonstrates the scattering cross 
section from an icosahedron enclosing a sphere. This is a complex, disjointed structure. 
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FIG. 18: (Color online) An application: SCS computation for a complex, multiply 
connected object (icosahedron enclosing a sphere), initially described using a point cloud. 

The structure is represented using a point cloud with Npt^ = 4888 points. Patches are 
constructed from this point cloud with varying g, to maintain an relative error of Er = 10~^ 
and basis functions of order p = 2 are constructed on the patches. The figure shows the 
surface current on the icosahedron and the SCS computed at = for an incident acoustic 
field along z. 

VII. CONCLUSION 



In this paper, we have developed and implemented a highly flexible framework for solving 
scattering from acoustically hard objects. The framework is very flexible in that the functions 
used to represent the fields are divorced from the underlying tessellation as continuity is built 
into representation space. This separation permits relatively easy modification of geometry 
and function representations, independently, so as to achieve convergence. Here, several 
benefits of the GMM has been demonstrated: namely, (i) ability to compute scattering from 
objects described using either a point cloud or a standard tessellation, (ii) ease of refining the 
patch size, order of surface, order of approximation and various combinations thereof, and 
(iii) ability to use mixtures of polynomial and non-polynomial functions. We are currently in 
the process of developing methods wherein this technique can be applied to solving Maxwell 
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equations, and the results will be presented elsewhere. 
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TABLE I: Patch -statistics for merged patches on parabohc surface 
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